Free energy of the Frohlich polaron in two and three dimensions. 
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We present a novel Path Integral Monte Carlo scheme to solve the Frohlich polaron model. At 
intermediate and strong electron-phonon coupling, the polaron self-trapping is properly taken into 
account at the level of an effective action obtained by a preaveraging procedure with a retarded 
' trial action. We compute the free energy at several couplings and temperatures in three and two 

dimensions. Our results show that the accuracy of the Feynman variational upper bound for the 
■ free energy is always better than 5% although the thermodynamics derived from it is not correct. 

Our estimates of the ground state energies demonstrate that the second cumulant correction to 
the variational upper bound predicts the self energy to better than 1% at intermediate and strong 
coupling. 
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I. INTRODUCTION 



The quasiparticle formed by an electron added to a polar crystal and the lattice deformation associated to it is 
called a "polaron" after Landau For non localized electrons in polar insulators and semiconductors, Frohlich et 
| al. ^introduced a model for a single electron interacting with dispersionless longitudinal optical phonons. Few years 
<—j . later, within the path integral formalism, Feynman proposed a very powerful variational procedure to compute the 
ground state energy and the effective mass of the model which provide accurate results at all couplings || . Extension 
i of the Feynman technique to compute properties at finite temperature appeared few years later by Osaka and 

more recently by Castrigiano et al. Q (see also Ref. ||). Recently the polaron problem has gained new interest as it 
C*~) . could play a role in explaining the properties of the high T c superconducting materials [Q. 

After a first Monte Carlo attempt by Becker et al. || , a true numerical test of the Feynman solution was provided 
by Alexandrou et al. [||. Applying an original extension of the Fourier Path Integral Monte Carlo (PIMC) method to 
' retarded actions, and special importance sampling technique, they provided exact (in the Monte Carlo sense) results 
, for the ground state energy and the effective mass of the Frohlich polaron in three dimensions. 

Recently the same model has been studied with a diagrammatic Monte Carlo (MC) technique. Attention has been 
focused on the ground state energy effective mass, phonon distribution and on the excitation spectrum at zero 
temperature [lljj . Also the effective mass for a lattice version of the Frohlich polaron in two dimensions has been 
obtained |||. 

In this letter we extend the numerical work on the Frohlich polaron to finite temperature and we provide results 
for the free energy in three and two dimensions. We focus on the intermediate and strong coupling regime which 
was only partially covered by the previous ground state MC investigations. Extrapolation of the free energy to zero 
temperature provides the ground state energy. The result of our analysis is that the Feynman variational method 
(FVM) at finite temperature provides estimates for the free energy within few percent at all couplings and dimensions, 
as it was already demonstrated for the ground state energy ||. However, due to its variational nature, FVM does not 
provide a consistent solution of the problem as thermodynamic derivatives of the free energy are not given correctly. 
In particular the entropy of the variational solution does not vanish with temperature as it should according to the 
third principle of thermodynamics. This is not just an academic issue since in the applications where the model could 
be invoked thermal effects are likely to be relevant. Indeed the phonon frequency in real materials is of the order of 
tenths of meV. 

We have formulated the problem in the framework of the discretized time PIMC method [ p"3[ . The high temperature 
density matrix has been obtained through an extension of the preaveraging technique (or cumulant method [|l3j) 
to retarded actions. This procedure regularizes the attractive coulomb-like divergence of the bare action at short 
distance (see below), and makes well defined and stable the discretized model fl5|| . Furthermore, in order to cope with 
the polaron self-trapping phenomena, we propose a suitable retarded trial action (RTA) around which the preaveraging 
procedure is more effective than the standard preaveraging with the free particle trial action (FPTA). At intermediate 
and strong coupling and at low temperature, the RTA improves the convergence with the number of time slices M 
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compared to the FPTA. Moreover, at given coupling, it removes a cross-over in the convergence of thermodynamic 
quantities with the number of slices and allows to obtain reliable extrapolations (M — > oo) from simulations with 
limited number of slices only. The methodological aspects of the present work are essential to attack the problem of 
the simulation of many interacting large polarons a system which has received some attention [l7[-|20|] as it could be 
relevant to describe the unscreened long range interactions in the CuO planes of the cuprates high T c superconductors 
pl[ . It is known that the free energy differences between different phases of the electron gas are quite small . The 
characterization of the electron-phonon effects on these phases requires an accurate method to compute the free energy 
for the many polaron system. The method presented here is similar to the one used by Alexandrou and Rosenfelder 
but in real rather than in Fourier space || . 

We start from the effective action of the Frohlich polaron, after the phonons have been eliminated with the usual 
path integral technique p3| . In polaronic units, Tiloto — 1, V^7 ( muJ Lo) — 1 (^lo is the phonon frequency and m 
the electron band mass) , and discretizing the imaginary time interval /3 in M intervals of width r it is 



S = S + S e - ph = - V / dhq 

2 l= l J(i-1)T 



(1) 



EE/,. dti ,. dh - 



2 V2^^7( 2 -i)r J(j-i)T |q(*0 - q(*2)| 

Here a is the coupling constant, q(i) is the electron position at imaginary time t and D(t,/3) = [e* + e^~ l \j\e^ — 1] 
is the phonon propagator. In the present work we consider the case of isotropic band mass (D = 3) and the case in 
which the band mass in one direction is infinitely large and the motion of the electron occurs on a plane (D — 2) 
The density matrix for the electron to go from r to r' in a time (3 is 

p(r,r'\/3) = J n^dq, J Dx^-^' e~ ST (2) 

where q., = q(jr), qo = q(0) = r, qnj = q(/3) = r', Xj(i) represent the paths of the j — th time slice (from q^-i to 
q^) and we have introduced a retarded trial action quadratic in the path coordinates 

St = So (3) 

^ M M r.i r 

+ o-EE^S / ^ / *2|q(il)-qfe)| 2 
1 l= lj=l J(i-l)r i(j-l)r 

where D™- = D(wr\i — j\,wf3) and C and w are variational parameters. In eq. (Q) (•■■) T . represents an average over 
paths belonging to the j — th time slice weighted by the trial action St- Using the Fourier representation of each time 
interval, and taking advantage of the gaussian character of the Fourier amplitudes which arises from the quadratic 
trial action (|^), we can analytically perform the preaveraging procedure |l4|,^5| to obtain an effective action between 
any pair of slices in terms of the end points of the slices only 



p(r,r'|/3)> Pc (r,r'|/3)= (4) 
e -<s-s T > Tz f £,„. e -s T 1 = 



N Md/2 . 
2-KT ) J 

For C = the standard free particle trial action local in imaginary time is implied and we get 

M 

So = r £|q«-q«-i| 2 (5) 

2=1 

8«* = -^ee r *• r * j( t; + v Bi) c (!£fLf\ ( 6) 
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where 



G(x) = erf{x) D = 3 (7) 

lxl (^)e-*'* D = 2 (8) 

= +u{q l -qi_i) - v(q,- - q,-_i) (9) 

£ 2 3 ( U ,«) = t[1-2( U 2 + i, 2 )] i^j (10) 

S 2 l ( U ,«) = 2r| U - W |(l-| U - U |) (11) 

In eq.(||), /of^) is the modified Bessel function p6[ ]. In the general case C / 0, the expression for £ e // is much 
longer and we will report it in a more extensive publication pq |. Here we just mention that its derivation is long but 
straightforward. 

We have sampled the path configurational space by Metropolis Monte Carlo method Jl3| with probability distribu- 
tion proportional to exp{— (So + S e ff)}. Even at strong coupling, the simple Levy reconstruction of paths |l3| was 
appropriate to efficiently move through configurational space. Typically between 20 and 400 time slices have been 
used in this study. The double integral over u and v in S e ff (see eq.(|6|)) is performed by Gauss method on the fly 

The method outlined above to derive the effective action is variational since F c = —j3 1 ln[Tr(p c )] > F = 
—f3~~ 1 ln[Tr(p)]/ (3 holds at any finite r, the equality being valid in the t — * limit. At given r, optimal trial ac- 
tion can be obtained minimizing the average of the effective action with respect to the parameters C and w in eq.(^). 
This is related to the excess free energy F^ x through the thermodynamic integration in the coupling constant a 

Fr(a) = U a da' <Seff / >a ' (12) 
P Jo a' 

We have checked in few significant cases that the optimal values for C and w are always close to the Feynman 
variational values B, w hich we have used in our simulations. 

We exploited eq.(|l2[) to obtain the free energy of the polaron. At given coupling and temperature, the coupling- 
constant integration has been performed by Gauss method. We have checked that order four in the Gauss integration 
is enough to obtain accurate results even at strong couplings. At each Gauss point a series of simulations for increasing 
number of slices were necessary to extrapolate the average effective action to r = 0. In figure |l| we show the typical 
behavior of the average effective action versus r at strong coupling and low temperature, and we compare data obtained 
with the FPTA (C = 0) and with the RTA (C ^ 0). The RTA improves the convergence at fixed r by a factor between 
two and three and it removes a crossover in r observed in the behavior for the FPTA. It allows to extract reliable 
extrapolations from data at large r only. This phenomenon can be understood comparing the diffusional properties 
of the paths with the two trial actions Pq] . The imaginary time mean square displacement A(t) =< |q(t) — q(0)| 2 > 
is the quantity appropriate to discuss the self-trapping of the polaron. Indeed, at strong coupling and at low enough 
temperature, the mean square displacement presents a crossover between a free particle-like regime at short time 
t < t x (A(t) = t((3 — t)/(3) and a localized regime at larger time t > t x (A(t) w constant) as it is shown in figure 
||. In the same figure we see that the results obtained with the RTA exhibit the self-trapping phenomena already for 
a limited number of slices (r > t x ) and they converge with M very quickly to the asymptotic behavior. Conversely, 
with the FPTA must be r < t x in order to correctly sample the self-trapped polaron. This explains the crossover in 
figure 0. 

In tables [| and [n| we report the values for the free energy at some selected thermodynamic points in three and two 
dimensions, respectively, and we compare with the results of the FVM both for the upper and lower bounds p7j| . The 
Feynman upper bound (FUB) is always close to the PIMC data (few percent). At intermediate coupling, FUB is more 
accurate in three rather than in two dimensions: the relative deviation of FUB from the PIMC data is maximum for 
a = 3 and D = 2 (between 2.5% and 4.5%) which corresponds approximately to a — 7 for D = 3 |24[ |. Data for D = 3 
and a — 5 suggest that at fixed coupling, the relative deviation is maximum around ftuho = 0.1 

At given a, the excess free energy for T < 0.1 is well reproduced by a polynomial function in terms of the even 
powers of T only. This is compatible with the third principle of thermodynamics as the entropy at T = must vanish. 
In figure ^ we show the excess free energy versus temperature at a = 9,0 = 10 and D — 3. We also report in the 
insets the entropy and the internal energy obtained as temperature derivatives of the polynomial function fitted to 
the PIMC data. Comparison with curves obtained from the FUB show the inadequacy of the variational method to 
provide correct thermodynamic derivatives of the free energy. 
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We estimate the ground state energy of the polaron by the coefficient of the order zero term in the fitting function 
for the excess free energy. The numerical values obtained for the considered cases are collected in table III . Our data 
compare well with a previous PIMC estimate H (Eq = —5.537(20) at a = 5, D = 3) and confirm that predictions 
based on the calculation of the second cumulant within FVM p8| are very accurate at any coupling. At D = 2 no 
such calculations has been performed. 

In conclusion we have developed an efficient scheme to simulate the single Frohlich polaron model within the 
discretized time Path Integrals Monte Carlo method. The new scheme is based on two key points: i) the extension 
of the preaveraging procedure to retarded actions; ii) the introduction of a suitable quadratic trial action, similar to 
the one used by Feynman, around which to perform the preaveraging. We have shown that the scheme is able to 
integrate the short (imaginary) time motion of the paths up to times larger than the self-trapping time of the polaron. 
Applying this new method we have obtained fully converged results for the free energy and the ground state energy 
of the polaron in two and three dimensions at intermediate and strong electron-phonon coupling. The present results 
show that the accuracy of the Feynman upper bound for the free energy is always better than 5%. 

Our work extends the one of Alexandrou and Rosenfelder Q which however formulated the problem in the framework 
of the Fourier Path Integral Monte Carlo. The advantage of the present formulation is that its extension to the 
many polarons system is straightforward. The present scheme can be used in connection with the recently proposed 
Restricted Path Integrals MC devised to circumvent the fermion sign problem Work in this direction is in 

progress. 



ACKNOWLEDGMENTS 



We acknowledge useful discussions and suggestions from DM. Ceperley, S. Fratini and J. Lorenzana. We acknowl- 
edge partial support from the MURST 1999 matching funds program. 



emails: sergio.ciuchi@aquila.infn.it, carlo.pierleoni@aquila.infn.it 
[1] L.D. Landau, Sov. Phys.3, 664 (1933). 

[2] H. Frohlich, H. Pelzer and S. Zienau, Philos. Mag. 41, 221 (1950), H. Frohlich, Adv. Phys. 3, 325 (1954). 
[3] R.P. Feynman, Phys. Rev. 97, 660 (1955). 

[4] Y. Osaka, Prog. Theor. Phys. 22, 437 (1959), J. Phys. Soc. Japan 21, 423 (1965). 

[5] D.P.L. Castrigiano and N. Kokiantonis, Phys. Lett. 96A, 55 (1983); D.P.L. Castrigiano, N. Kokiantonis and H. Stierstorfer, 

Phys. Lett. 104A, 364 (1984). 
[6] D.C. Khandekar and S.V. Lawande, Phys. Rep. 137, 115 (1986). 

[7] "Materials & Mechanism of Superconductivity, High Temperature Superconductors V" , Y-Sheng He, Pei-Heng Wu, Li- Fang 

Xu and Zong-Xian Zhao eds, Physica C 282-287 (1997). 
[8] W. Becker, B. Gerlach and H. Schliffke, Phys. Rev. B 28, 5735 (1983). 

[9] C. Alexandrou, W. Fleisher and R. Rosenfelder, Phys. Rev. Letts 65, 2615 (1990); C. Alexandrou and R. Rosenfelder, 
Phys. Rep. 215, 1 (1992). 
[10] N.V. Prokof'ev and B.V. Svistunov, Phys. Rev. Lett. 81, 2514 (1998). 

[11] A.S. Mishchenko N.V. Prokof'ev, A. Sakamoto and B.V. Svistunov, Phys. Rev. B 62, 6317 (2000). 
[12] P.E. Kornilovitch, Phys. Rev. Lett. 81, 5382 (1998). 
[13] D. M. Ceperley, Rev. Mod. Phys. 67, 279 (1995). 

[14] J. D. Doll, D. L. Freeman and T.L. Beck, Adv. Chem. Phys. 78, 61 (1990). 
[15] S.Ciuchi, JLorenzana and C.Pierleoni, Phys. Rev. B 62, 4426 (2000). 

[16] D.M. Ceperley, "Path Integral Monte Carlo methods for fermions", in Monte Carlo and Molecular Dynamics of Condensed 

Matter Systems, Ed. K. Binder and G. Ciccotti, Editrice Compositori, Bologna, Italy, 1996. 
[17] L.F. Lemmens, J.T. Devreese and F.Brosen, Phys. Stat. Sol. B 82, 439 (1977). 
[18] G. D. Mahan, Many Particle Physics (Plenum, New York, 1990). 
[19] G. De Filippis, V. Cataudella and G. Iadonisi, Eur. Phys. J. B 8, 339 (1999). 
[20] S. Fratini and P. Quemerais, Eur. Phys. J. B 14, 99 (2000). 

[21] S. Lupi, P. Maselli, M. Capizzi, P. Calvani, P Giura and P. Roy, Phys. Rev. Letts. 83, 4852 (1999). 
[22] M.D. Jones and D.M. Ceperley, Phys. Rev. Letts. 76, 4572 (1996). 
[23] R. P. Feynman, Statistical Mechanics, Benjamin (New York, 1972). 



4 



[24] Wu Xiaoguang, F.M. Peeters and J.T. Devreese, Phys. Rev B 31, 3420 (1985). 
[25] J.T. Titantah, S. Ciuchi and C. Pierleoni, to be published. 

[26] M. Abramowitz and I. Stegun, Handbook of Mathematical Functions, Dover Publications (New York, 1972). 

[27] Within Feynman method a lower bound for the free energy can be obtained as Ft + fi^ 1 < S — St >< F where Ft is the 
free energy of the quadratic Feynman action and the average is performed with the true polaron action S. Since F T and 
— St are both negative quantities the lower bound is obtained when Ft = St = and it is simply < S > /f3 independent 
on the trial action. 

[28] J.T. Marshall and L.R. Mills, Phys. Rev B 2 3143 (1970). 



if) 
< 



0.1 



0.01 









tx 


• 







0.01 



0.1 



1 



FIG. 1. a = 9, P = 10. Relative deviation of the excess action with respect to its extrapolation at r = 0, 
ASeff = [< S e ff >t — < S e ff >o]/ < S e ff >o, versus t. Comparison between FPTA data (squares) and RTA data 
(circles). Lines are power law fits to the data. The crossover t x discussed in the text is clearly seen in the FPTA data. 
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FIG. 2. a = 9, (3 = 10. Open (closed) symbols are data obtained with the FPTA (RTA). Data for M = 16 (squares) and 
M — 32 (circles) are reported. The thick continuous line is the results of a fully converged simulation (M = 256 with the RTA), 
while the thin continuous line is the prediction of FVM. The dashed line is the free particle behavior while the dashed vertical 
line indicates t x . 
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FIG. 3. Excess free energy versus temperature at a = 9. Closed circles are PIMC data and the continuous line is a polynomial 
fit to the data. The dashed line is the FUB. In the inset A we show the entropy vs temperature as obtained by deriving with T 
both the polynomial fit to the PIMC data (continuous line) and the FUB (dashed line) . In the inset B we report the internal 
energy F + TS. 



TABLE I. 


D=3. Comparison between the PIMC free energy and the upper bound (FUB) and the lower bound (FLB) of the 


Feynman variational method. 








a 


13 


FLB 


PIMC 


FUB 


1.0 


1.0 


-1.87001 


-1.838(2) 


-1.830 


1.0 


5.0 




-1.103(1) 


-1.095 


1.0 


10.0 


-1.10077 


-1.053(1) 


-1.043 


3.0 


1.0 


-5.90912 


-5.659(1) 


-5.638 


3.0 


5.0 


-3.70067 


-3.479(2) 


-3.419 


3.0 


10.0 




-3.251(3) 


-3.244 


5.0 


1.0 


-10.8646 


-9.827(5) 


-9.696 


5.0 


5.0 


-7.16466 


-6.12(1) 


-6.018 


5.0 


10.0 


-6.71199 


-5.800(5) 


-5.678 


5.0 


15.0 


-6.49491 


-5.654(3) 


-5.590 


5.0 


20.0 


-6.41828 


-5.593(2) 


-5.550 


5.0 


40.0 


-6.37900 


-5.539(2) 


-5.493 


9.0 


10.0 


-17.94248 


-12.43(1) 


-12.13 


9.0 


15.0 


-17.37586 


-11.94(1) 


-11.91 


9.0 


20.0 


-17.33489 


-11.85(1) 


-11.80 


9.0 


40.0 


-17.32615 


-11.80(2) 


-11.64 


TABLE II 


D=2. Comparison between the PIMC free energy and the upper bound (FUB) and the lower bound (FLB) of the 


Feynman variational method. 








a 


P 


FLB 


PIMC 


FUB 


1.0 


10.0 


-1.79381 


-1.6940(6) 


-1.6778 


1.0 


15.0 


-1.76034 


-1.6663(9) 


-1.6576 


1.0 


20.0 


-1.74893 


-1.6535(7) 


-1.6483 


1.0 


40.0 




-1.6412(9) 


-1.6354 


3.0 


10.0 


-7.66997 


-6.044(5) 


-5.777 


3.0 


15.0 


-7.45560 


-5.771(6) 


-5.671 


3.0 


20.0 


-7.40605 


-5.716(3) 


-5.620 


3.0 


40.0 




-5.677(6) 


-5.547 
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TABLE III. Polaron ground state energy. Comparison between our present estimates (PIMC), variational upper and lower 
bounds (FUP and FLB), data from ref. [27] for D = 3 and from ref. [23] for D = 2 (PT). 



D 


a 


FLB 


PIMC 


FUB 


PT 


3 


5.0 


-6.360 


-5.523(3) 


-5.440 


-5.52 


3 


9.0 


-17.292 


-11.794(2) 


-11.49 


-11.7 


2 


1.0 


-1.735 


-1.635(2) 


-1.616 


-1.63477 


2 


3.0 


-7.402 


-5.670(3) 


-5.467 


-5.28812 
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